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ABSTRACT 


The  one  dimensional  problem  of  the  diurnal  variations  of  moisture  and 
temperature  at  different  depths  of  the  thermally  active,  non-frosen, 
non-swell ing  soi 1  layer  is  considered. 

It  is  assumied  that  we  knowi  (ij  values  of  all  meteorological  elements  at 


anemometer  level  at  any  time^  f2J  Moisture  and  temperature  distribution  with 


depths  at  some  initial  moment.  ^3)  Physical  parameters  and  caracter i st i cs  of 


the  soil  and  vegetation  cover  which  we  shall  call  canopy  for  brevity. 

The  model  consists  of i  (ij  Partial  differential  equations  of  the 
subsurface  hydrology  for  a  non-saturated  vapor-liquid  flow  in  the  soil,  allowing 
for  the  moisture  removal  by  the  roots  of  plants.  (2)  Algebraic  empirical  and 
ordinary  differential  equations  suggested  by — Drar dorff  <197Q^for  the 
description  of  the  radiative  and  turbulent/heat  transfer  between  the  atmosphere 


and  the  soil  through  the  canopy  allowing  for  the  influence  of  condensation  or 
evaporation  at  the  foliage  and  also  the  retention  of  a  part  of  the  precipitation 
by  its  the  foliage,  (s)  Moisture  and  vapor  balance  equations  at  the  surface  of 
the  soil.  (4)  Continuity  conditions  for  heat  and  vapor  fluxes  at  the  top  of  the 
canopy,  (s)  Semi emp i r i cal  algebraic  equations  of  the  type  of  surface  layer 
equations  for  a  thin  atmospheric  layer  between  the  canopy  and  anemometer  level. 

The  method  of  solution  is  the  following: 

1.  Transforming  and  solving  the  equation  describing  the  effect  of  canopy,  we 
construct  boundary  condition  for  moisture  and  temperature  at  the  soil  surface 
which  takes  into  account  that  effect.  2.  The  soil  temperature  and  moisture 
equations  as  well  as  the  boundary  and  initial  conditions  are  approximated  by  the 
finite-difference  equations  in  accordance  with  known  Crank-Nicol son  schemce 


respectively.  These  schemes  are  quite  simple  and  provide  the  convergence  and 
the  second  order  ot  the  accuracy  with  respect  to  the  time  and  space  steps.  3. 


The  -f  i  n  i  te-di  t-f  erence  equations  are  solved  by  the  tridiagonal  algorithm  in 

I,  p.  I 

combination  with  iteration  method.  ' 

The  solution  has  been  programmed.  Resuites  of  computations  are  in 


agreement  with  experimental  data.  A  few  examples  of  soil  moisture  computations 
are  presented  together  with  experimental  data. 
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1.  INTRODUCTION 

Information  on  soil  moisture  and  temperature  at  the  surface  and  at 
different  depths  is  of  interest  in  numerous  fields  of  human  activity,  e.g., 
underground  structures  (water  pipes,  sewage,  electric  and  other  cables),  road 
building,  some  military  problems  (for  example,  road  practicability  for  tanks), 
etc.  In  agriculture  and  agrometeorology,  this  information  is  often  very 
necessity.  For  semiarid  and  arid  zones  such  information  becomes  especially 
important . 

There  is  usually  a  great  deal  of  technical  difficulty  and  expense  involved 
in  conducting  special  soil  moisture  and  temperature  measurements.  Therefore,  it 
is  very  tempting  to  have  a  model  which  would  allow  the  computation  of  soil 
moisture  and  temperature  from  routine  meteorological  observations  available  for 
many  regions  of  the  globe.  Moreover,  such  a  model  would  enable  the  prediction 
of  soil  moisture  and  temperature  using  the  meteorological  forecast. 

The  construction  of  a  model  encounters  some  difficulties  such  as 
inhomogeneities  of  the  land  surface  and  of  the  soil  at  various  depths,  lack  of 
information  on  physical  parameters  and  presence  of  different  urban  constructions 
and  artificial  covers  (such  as  concrete),  as  well  as  natural  vegetation  cover. 

It  has  been  shown  that  vegetation  can  significantly  influence  the  heat  and 
moisture  transfer  between  the  earth's  surface  and  the  atmosphere.  Therefore, 
since  most  of  the  land  surface  of  the  globe  is  covered  by  vegetation,  it  is 
desirable  to  include  the  vegetation  effect  in  a  model  which  simulates  the  heat 
and  moisture  transfer  at  the  land  surface. 

The  object  of  the  present  paper  is  to  construct  a  model  for  an  operational 
soil  moisture  and  temperature  computation  from  meteorol ogi cal  data,  taking  into 


account  the  et-fect  o-f  the  vegetation  cower  which  hereatter  will  be  re-ferred  to 
as' a  canopy.  To  avoid  the  other  di -f-f  i  cul  t  i  es  mentioned  above,  we  assume  that 
the  soil  is  homogeneous,  and  the  land  surface  is  horizontal  and  -free  o-f  urban 
constructions. 

To  our  knowledge,  there  are  no  publications  devoted  even  to  such  a 
simplified  problem,  although  there  are  many  papers  in  which  soil  moisture  and/or 
soil  temperature  computations  have  been  carried  out.  These  works  have  been 
directed  at  the  investigation  of  soil  physics  or  of  a  methodology  of  numerical 
solution  <e.g.,  [1],  [2]).  In  these  papers,  the  boundary  conditions  at  the 
surface  are  mostly  simplified,  e.g.,  evaporation  from  the  surface  is  assumed 
constant  and  known.  Such  assumptions  are  not  sufficient  for  our  problem. 
Included  in  this  group  are  some  papers  in  which  an  attempt  is  made  of  coupling 
the  soil  part  of  the  problem  with  a  rather  primitive  atmospheric  model  assuming 
bare  soil  with  attention  directed  mainly  at  the  soil  <e.g.,  [3],  [4]>.  These 
works  have  allowed  us  to  arrive  at  the  important  conclusion  that,  in  principle, 
soil  moisture  and  temperature  computat i on^ wi th  the  accuracy  necessary  for  many 
aims, i s  possible . 

Another  group  of  papers  is  devoted  to  meteorological  and  cl imatolooicil 
problems  in  which  the  soil  moisture  and  temperature  play  an  auxiliary  role.  In 
the  majority  of  these  works  the  soil  part  of  the  problem  is  highly  simplified 
(e.g.,  [5])  or,  if  not  oversimplified,  cumbersome  calculations  and  the  lack  of 
necessary  Information  make  It  very  difficult  to  reliably  estimate  the  accuracy 
with  which  soil  moisture  is  determined  (e.g.,  161).  In  this  group  of  papers,  we 
refer  only  to  those  in  which  attempts  to  allow  for  the  vegetational  cover  effect 
have  been  made.  Obviously,  canopy  models  accessible  for  inclusion  into  complex 


models  containing  soil  and  atmospheric  parts,  must  be  sufficiently  simple  with 


respect  to  mathematics.  Nevertheless,  there  have  been  attempts  to  construct 
rather  complicated  two-level,  or  even  multi-level,  biosphere  models  tor 
inclusion  into  general  circulation  models  ot  the  atmosphere  [71.  Deardortt 
<1978)  [5]  has  suggested  a  single-level  canopy  model,  sut t i c i entl y  simple  with 
respect  to  mathematics  and  at  the  same  time  rather  comprehensive  with  respect  to 
physics.  The  model  is  based  on  moisture  and  heat  balance  equations  tor  the 
vegetation  layer  and  tor  the  surtace  ot  the  earth  and  on  some  empirical 
relationships.  Deardortt's  model  has  already  been  employed  by  some  authors 
<[<$],  [8]).  In  the  present  paper  we  also  use  Deardortt  s  model  while 
introducing  minor  modi t i cat i ons  tor  its  improvement. 

Unlike  Deardortt  and  the  authors  mentioned  above,  we  have  made  some 
preliminary  analytical  transtormat i ons  in  an  attempt  to  simplity  the  solution  ot 
the  vegetational  part  ot  the  problem.  The  present  paper  includes  only  the 
physical  and  mathematical  tormulation  ot  the  problem  and  its  method  ot  solution. 
Results  ot  computations  and  comparisons  with  observations,  as  well  as  physical 
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2.  NOTATION 

Arab  ic 

a  anemometer  level  height 

by  empirical  value  dependent  on  type  (textural  class)  of  soil  only 

C  volumetric  heat  capacity  of  soil 

C|  heat  or  moisture  transfer  coefficient  for  the  foliage  element 

heat  and  moisture  transfer  coefficient  applicable  to  bare  soil. 
C||^!  to  soil  under  a  canopy,  to  the  top  of  a  dense  canopy, 

moisture  diffusivity  of  soil 
Dy  thermal  moisture  diffusivity  of  soil 

E  evaporation  rate,  if  E>0,  or  condensation  rate  is  E<0 

E^^  transpiration  rate  (per  unit  soil  surface  area) 

g  gravity  acceleration 

h  thickness  of  the  vegetation  layer 

H  sensible  heat  flux,  positive  upwards 

i,J  depth  and  time  indices 

K  hydraulic  conductivity  of  soil 

L  latent  heat  of  vaporization 

Ly,  length  of  roots  per  unit  soil  volume 

N  number  of  levels  in  soil 

P  precipitation  rate 

Q  soil  heat  flux  density,  positive  downward 

q  spec i f i c  humi di ty 

‘’5 


r 


saturated  humi di ty 

fraction  of  potential  evaporation  rate  from  foliage 


'vw 


t 

T 

u 

U 

z 

2. 


Greek 

cl. 


A  z,  At 


'«l<w 


,9- 


atmospheric  resistance 

resistance  coe-f-f  i c  i ent  dependent  upon  plant  type 
generalized  stomatal  resistance 
sur-face  runo-f-f 

downward  directed  longwaue  radiative  'flux 

downward  directed  shortwave  radiative  -flux 

maximum  value  of  Rjj 

gas  constant  for  water  vapor 

bulk  Richardson  number 

t  ime 

temperature 
wind  speed 

flux  density  of  soil  water  transfer,  positive  downwards 
vertical  coordinate  positive  downwards  (depth  in  soil) 
lower  boundary  of  the  solution  domain 
depth  of  root  zone 

fol i age  albedo 
ground  surface  albedo 
depth  and  time  increments 
volumetrit.  moisture  content 
residual  water  content 

mass  of  liquid  water  retained  by  foliage  per  unit  horizontal  area 
maximum  value  of  ,  beyond  which  runoff  to  soil  occurs 

minimum  value  of  "O'  in  the  root  zone 


7 


value  o-f  V  at  saturation 
i/iilting  point  value  of 
dens i ty  of  water 
dens i ty  of  air 

thermal  conductivity  of  soil 
soil  moisture  potential 
Stefan-Bol tzman  constant 
cloud  fraction 

foliage  shielding  factor  of  ground  from  shortwave  radiation 
(area  average) 
fol i age  emi ss i v i ty 
ground  surface  emissivity 

Subscr ip  ts 

a  reference  "anemometer  level"  height 

f  f ol i age  surface 

g  value  at  the  ground  surface 

af  mean  value  within  a  canopy 

h  value  at  the  top  of  the  canopy 

s  value  for  saturated  soil 


3.  FORMULATION  OF  THE  PROBLEM 


We  consider  the  problem  o-f  diurnal  moisture  and  temperature  variations  in 
the  upper  thermally  active  soil  layer  which  we  assume  to  be  non-swelling.  In 
addition,  we  assume  that  we  know  the  values  o-f  the  wind,  temperature,  pressure, 
humidity,  precipitation,  direct  and  scattered  radiation  at  the  shelter  level  at 
all  times  and  that  we  have  all  the  necessary  in-formation  about  the  physical 
parameters  of  the  soil  and  the  canopy.  We  imply  that  there  is  a  layer  o-f  air 
between  the  shelter  and  the  canopy. 

To  describe  the  heat  and  moisture  trans-fer  in  the  soil,  we  employ  a  system 
of  subsurface  hydrology  equations  for  unsaturated  vapor-liquid  flow  191,  IlOI. 

We  assume  that  nothing  depends  on  the  horizontal  coordinates  and  that  the  ground 
surface  is  horizontal.  The  soil  moisture  equations  will  be  supplemented  by  the 
water  extraction  term  till  which  parameterizes  the  moisture  removal  from  the 
soil  by  the  roots  of  the  plants. 

As  a  result,  the  soi 1 -moi sture  equations  are  reduced  to  the  following 


one-dimensional,  non-stat i onary  system  of  partial  differential  equations. 


i 


i 


In  these  equations  the  z  axis  and  -fluxes  are  positive  downward. 


Unlike  [9],  eq.  <4)  does  not  take  into  account  the  trans-fer  o-f  latent  heat 
by  vapor  movement  induced  by  the  moisture  gradient.  The  possibility  o-f  omitting 
this  et-fect  can  be  shoixn  by  elementary  estimation.  The  -fact  is  that  the  thermal 
vapor  di -f-fusi  V  i  ty  coefficient  becomes  extremely  small  for  soil  moisture  values 
observed  in  nature  <see  e.g.,  Fig.  2  from  [9]). 

Neglecting  the  effect  of  hysteresis,  one  can  consider  parameters  C,  D^, 

Qj.  >  •<>  ^  and  also^  as  unique  functions  of  and  of  the  textural  class  <or 
type)  of  the  soil,  whereas  parameters  ^  i  ^  depend  only  on  textural  class 
of  the  soil.  Following  [123-C15]  we  set^ 

+0  063(-'l'y  ( axje/sec-Cft,-°c) 

The  values  ^  I  and  S-j- for  different  soil  are  given  in 

Tabl e  2  from  112]. 

To  parameterize  the  soil  moisture  removal  by  roots,  we  adopted  an 
expression  for  S,  suggested  recently  by  [11]  (see  also  [161). 


V  The  expression  for  A  approx imates  quite  accurately  the  experimental 
data  plotted  in  Fig.  4  from  [131. 


S 


)KE^/  («> 

assuming  that  and  z^  are  known. 

We  now  proceed  to  the  formulation  of  boundary  and  initial  conditions  of  the 
problem. 

Since  the  ground  surface  is  assumed  to  be  horizontal,  one  can,  with  slight 
simplification,  consider  that  surface  runoff  occurs  only  if  the  precipitation 
water  has  no  time  to  infiltrate  into  the  soil  and  to  evaporate  from  the  ground 
surface,  ighich  becomes  saturated  instantly: 


either  and  R  >  0 

or  ^  0^^  and  R  =  0 


at  z  =  0 


where 


(7) 


(8) 


(9) 


(7)  is  actually  an  upper  boundary  condition  for  the  soil  moisture.  <8)  is  the 
surface  moisture  balance  equation.  <9)  allows  for  the  retention  of  part  of  the 
precipitation  by  the  foliage.  As  an  upper  boundary  condition  for  the  soil 
temperature,  we  shall  use  the  ground  surface  heat  balance  equation 


-  11 


where 


-^c*iS.i(i-'S-,)(^^^‘]^T^ 

C an  -o.i'fS /a,  B-^<os9, 


o.f‘/ 


9.  -s  a.s-^ 


Ue  assume  that  R-,^  s  known  -for  any  time  -from  observations.  It  is  also  possible 
to  use  astronomic  formulas  tor  ,  as  in  t63. 

Eq .  (10)  contains  additional  terms,  which  take  into  account  the  ettect  o-f 
the  canopy.  It  0,  eq.  (10)  becomes  the  well-known  ground  surtace  heat 

balance  equation  tor  the  bare  soil. 

Equations  which  have  to  parameterize  the  intluence  ot  the  canopy  and 
meteorological  conditions  on  the  heat  and  moisture  transter  between  the 
atmosphere  and  the  ground  are: 

1.  Canopy  energy  balance  equation 


('Sr 


12 


where 

“  j2  ‘^■Vj  ^<1/  (%  ■  %f ) 


C«  =  7 1-10  ^{‘U.^^+  03  H,/sn ) 


2.  The  relationship  between  the  moisture  characteristics  at  the 
and  the  in-canopy  air  humidity. 


where 


=i-s.4^4^  A’l^) 


dfvfci^ 


/ 


<1  3) 


(14) 


-fol  i  age  sur-face 


<15) 


<14) 


/  OJf  there  is  condensation  onto  the  -foliage  sur-face 

9.!^  %  CD 

S'cW 

I  1  it  there  is  evaporation  -from  the  -foliage  sur-face 


^  r  ^Sr>7ccX _ 


( 


(^7-j  ^■exp(i7.%6^ 


T-3r.i€  J 


This  empirical  tormula  is  actually  one  possible  version  o-f  the  well-known 
Cl ausi us-Cl apeyron  equation. 

3.  The  conservation  equation  -for  : 


< 


where 


4.  Continuity  conditions  ot  the  heat  and  water  vapor  tluxe*  at  the  top  ot  th 
canopy  layer 


<22) 


where 


<23) 


We  denote 


<24) 


The  last  two  relationships  are  constructed  in  accordance  with  an  idea  ot  the 
simplest  linear  interpolation  wi/A  Cj  ,  also  suggested  in  t5]. 

5.  Continuity  condition  for  q  at  the  ground  surface,  g. 


<25) 


6.  Semiempirical  relationship  for  C,.  (see  [0]) 

Ht, 


Cu  = 


^.X-lo^yO -f  US R;^]  ^  cf 


where 


2/  ^  Cf  ^  ^ 

i.  i  ni/^ec  if  S, 


(26) 


15 


These  relationships  imply  that  a  thin  constant  -flux  layer  exists  betiMeen  the 
canopy  and  the  anemometer  level  . 

For  the  formulation  of  the  lower  boundary  conditions  for^  and  T  ,  we 
assume  that  at  a  certain  depth  (which  is  usually  1-2  m) ,  diurnal  variations  of  ^ 
and  can  be  neglected.  We  consider  two  casesi  the  soil  below  that  depth  is 

dry;  the  soil  below  that  depth  is  saturated  (water  table): 

Thus,  the  loiwer  boundary  conditions  for  ^  and"7*i 
•0*  ^  ft  =  const  (dry  soil) 


are 


or 


=  const  (saturated  soil) 


at  X-Z, 


(27) 


7”  ~  const 

The  values  1  and  are  assumed  to  be  known  from  climatic  and  edaphic 

data. 

To  close  the  problem,  one  needs  to  formulate  the  initial  conditions. 

We  assume  that,  at  a  given  moment  of  time  't-  0,  ^  and  T  are  given 
funct i ons  of  X  i 


&  =  d-,Cs)  ,  -r^Tjz)  fo<z<z) 

^  ^  (28) 

Eq.  (20)  also  requires  an  initial  condition  for  ^  .  Evidently  it  is 
difficult  to  obtain  information  concerning  liquid  water  on  the  foliage. 
Therefore  we  set 

^  =^0  at  "d  =  0  (29) 

assuming  that  our  computation  starts  at  the  moment  when  all  the  water  retained 
by  the  foliage  after  the  last  rain  has  evaporated. 


4.  SOLUTIOt^  OF  THE  PROBLEM 


First  o-f  all  we  introduce  a  ^  i  n  i  te-di -f  •ference  time  gridlL®jA^  , 

J  J 

J=  1,2 . where  tp  is  any  function  or  parameter  of  the  problem  (including 

those  which  are  independent  of  Z  ) . 

To  clarify  the  explanation  of  the  computation  scheme,  we  shall  assume  that 
are  known.  We  then  describe  the  iterations  involved  in  advancing  to  the 

J'* 

next  time  step.  This  procedure  will  then  be  used  to  continue  for  all  j  .  For 
brevity  we  shall  omit  the  index  when  it  isj  .  This  notation  permits  us  to 
regard  eqs.  (12)-(19),  <21)-<24),  as  written  for  the  moment  . 

Eq .  (20)  can  be  approximated  by  the  following  finite-difference  form 

This  implicit  scheme  is  chosen  in  order  to  refer  dependent  variables  in 
(30)  to  the  moment  'tj  ,  as  in  eqs.  (12)-(19),  (21)-(24). 

Considering  the  crudeness  of  eq.  (20),  the  first  order  approximation  of  the 


time  derivative  in  (30)  appears  to  be  sufficient. 


Let  us  assume  temporarily  that  and 


^  (i  .e.,  And 


)  are  knwjn, 


Then,  as  one  can  show,  (12)-(19),  (21)-(24)  together  with  (30)  is  a  closed 
algebraic  system  from  which  we  shall  find  in  the  following  way. 

First,  substituting  (13),  (23)  and  (24)  into  (22)  and  making  some  simple 
transformations,  we  obtain 


13 


c^/[c,  Kf  ^  ^ 


In  the  right  hand  sides  of  <3D  and  (32)  all  quantities  with  the  exception 

o-f  ^  and  T  ,  are  known.  Indeed  O  and  7”  are  known  -from  obser  uat  i  ons , 

y  i  /A  A  '  ^ 

and  ^  are  known,  as  a  consequence  of  our  assumption  and  eq.  (25)  respectively. 

Concern i ng  and  ,  we  note  the  following.  In  [5]  it  was  assumed 

that  is  known  constant.  We  calculate  from  (26) ,  but  in  (  )  we 

shall  use(^"7^J‘.  .  Then  and  6? .  can  be  cal  cul  ated  wi  th  the  aid  of 

7  *  f  7 

(33).  It  should  be  noted  that  (31)  and  (32)  are  similar  to  relationships 

suggested  in  [51,  where  C(^  ,  0(.f  and  were  assumed  to  be  given  functions  o-f 

*  f 

alone.  Eqs.  (15)  and  (31)  can  be  solved  for  and  .  Substituting  ^ 

obtained  in  this  way,  into  the  expression  for  from  (13),  we  get 


,«<here 


^  '  fj]’’ 


f-lote,  that  O^Ti^l  andO^^,  ^1.  Therefore  the  sign  of  £4  depends  only 


the  sign  the  di-f-ference  in  the  brackets  in  (34). 


Substituting  (34)  into  (12)  yields  an  equation,  which  links'^  with/*  ! 


where 


in  which 


7^^  <3?> 

The  expressions  i.T\6  CAn  be  treated  as  given  -functions  of  their 


arguments. 


From  (34)  and  (36)  we  obtain 


Ef 


In  this  formula,  unlike  (34),  the  explicit  dependence  /y on is  absent. 


Combining  (13),  <31),  (35)  and  (40),  we  derivei 


Recalling  the  expressions  for  :  (13),  (34)  and  (40),  one  can  understand  from 
(41)  why  the  differences  q^<'^)-q,  q^'^)-q^and  q^cy-q^^always  have  the  same  signs. 
Thus,  from  (34)  and  (36)  it  follows  that  the  sign  of  to  de  term i nes 


9 


whether  condensation  or  evaporation  is  occur  ing  at  the  -foliage: 


H  £(ta<o  ~~  condensat  i 


H  -  evaporatii 


account  of  <1<S)  and  <17)  in  the  case  of  condensation,  can  be  found 


from  eq.  (36)  in  which  r  =  1. 


According  to  (17),  and  (21)  =  0  in  this  case.  Therefore,  in  the  case 


of  condensation  (30)  is  reduced  to 


In  the  case  of  evaporation  from  the  foliage,  an  equation  for  can  be 
obtained  by  excluding  r  and  ^  from  (16),  (30)  and  (36). 

Indeed,  in  this  case,  solving  (16)  and  (36)  for  ^  and  r*  respectively,  and 
allowing  for  (41),  we  express  and  K  in  terms  of 


Substituting  (40),  (44)  and  (45)  into  (21)  yields  an  expression  for  the 


transpiration  rate 


=  s;  a.^  \fa  f xl -Cf  l/ri  -  Of) 


(46) 


which  accompanies  evaporation  trom  the  -foliage. 


If  we  substitute  <40),  <44)  and  (4<S)  into  <30),  we  arrive  at  an  algebraic 
transcendental  equation  for  ~T^  in  the  evaporation  case 

where  the  functions  ^  and  F  are  specified  by  <44)  and  (36)  respectively  and 


efnycL^ 


Eq.  (36)  at  1”  =  1,  as  well  as  at  T  y  solved  numerically  by  the 

riewton-Raphson  method  C17].  Eq .  (47)  is  solved  using  the  so-called  "rule  of 
false  pos i t i on"  C 17] . 

As  a  result  of  (42),  (44),  (45)  and  (47)  one  can  find  and  determine 
whether  condensation  or  evaporation  at  the  foliage  is  occurring  at  the  moment'^- 
Details  are  given  in  APPENDIX  A. 

Once  we  obtain  ,  then,  using  (13),  (32),  (37),  (40),  (41),  (43)-(46),  we 
can  calculate  the  quantities: 


,r  q  q  T  E,  and  , 


which  are  necessary  for  the  computation  of  the  soil  temperature  T(z,t)  and 
mo i stur e  ^(z.t)  . 


I'l  «■— »  1'."  I'l.'  M"]  rrf’Ji  Wg'i  ti-w.  ■  'j 


-2i- 

E1  imi  nat  i  ng  IJ,  jj  and  GL -f  r  om(  1  )-(4)  yields 

ciZ  .  i 

7)1  "02 


';50) 


&3Z 


'dZ 


-»z^r  -dx 


<51) 


Note  that  these  equations  are  written  in  their  "prior  to  introducing  the 
time  indices"  form. 

Let  us  consider  the  soil  temperature  eq.  (50).  Sophocleous  (1979)  [181 
solved  the  temperature  equation  in  a  more  general  torm  (including  the  moisture 
in-fluence  term).  His  estimates  indicated  that  the  dependence  o-f  the  soil 
temperature  on  the  soil  moisture  tield  is  weak. 

This  tact  justified  the  use  of  the  simplified  eq.  (4)  and  also  shows  that 
changes  of  the  parameters  C  and  A  with  "9^  only  slightly  affect  T.  This  is 
reasonable  because  the  changes  of  C  and  A  have  to  be  mutually  compensated  to  a 
certain  extent. 

The  variables  (49)  depend  on  ft  only  through  the  soil  moisture  potential 

a 

Aj)  ,  which  appears  in  (25).  It  is  easy  to  show  that  the  exponential  factor  in 
<25>  is  close  to  unity,  if  the  soil  is  even  slightly  wet. 

However,  in  the  case  of  extremely  dry  soil,  the  value  of  ft  is  very  small. 

9 

Therefore,  in  both  cases  the  effect  of  ^  on  the  values  (49)  is  insignificant. 

Thus,  taking  into  account  the  weak  dependence  of  the  soil  temperature  field 
on  the  soil  moisture,  we  shall  use  j  instead  of  C^Jj  solving  eq. 

( 50 ) . 

To  transform  (10)  into  a  convenient  linear  boundary  condition  for  T,  we 


shall  use  the  fact  that  the  numerical  procedure  for  solving  eq.  (50>  contains  a 


certain  iteration  process. 

Let  us  turn  to  an  approximate  relationship 


';52) 


which  is  valid  when  the  temperature  difference  between  two  consecutive 
iterations  is  relatively  small. 


In  ( 52)  is  taken  from  the  (n  -  l)th  iteration.  In  analogy  with  the 


time  indices,  hereafter  the  iterative  index  in  brackets  will  be  omitted  vjhenever 


it  is  (n)  (i.e.,  referring  to  the  n-th  iteration).  For  example  in  (52) 


/t 

means  (T.).  . 


} 


Substituting  (13)  and  (52)  into  (10),  we  arrive  at  a  linear  boundary 
condition  for  T  at  the  ground  surface. 


where 


(53) 


(54) 


Eq .  (50)  is  rewritten  in  the  finite  differences  in  accordance  with  the 
we  11  known  Crank-Mi  col  son  scheme.  The  boundary  and  initial  conditions  ((27). 
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<53)  and  <28))  are  also  presented  in  -finite  difference  form.  Such  a 
finite-difference  problem  is  solved  at  each  iteration  with  the  aid  of  the 
so-called  tridiagonal  algorithm  (see  e.g.  [19]). 

Details  are  presented  in  APPENDIX  B. 

Once  T(z,t)  is  known,  we  proceed  to  the  solution  of  the  soil  moisture 
equation  (51)  with  the  boundary  and  initial  conditions  (7),  (27)  and  (28). 

For  this  we  make  use  of  the  predictoi — corrector  equations  (suggested  by 
Douglas  and  Jones  (19<43)  [20]  and  recommended  in  [19]  for  the  type  of  problems 
under  consideration)  again  solving  them  with  the  aid  of  tr i di zgnnal  algorithm  at 
each  iteration.  By  including  a  special  logical  scheme  in  the  iteration  process 
(see  Fig.  4),  we  provide  the  satisfaction  of  the  surface  boundary  conditions 
(7) . 

Details  are  presented  in  APPENDIX  C, 

Fig.  1  shows  a  simplified  flowchart  of  the  soil  moisture  and  temperature 

calculation  at  given  time  step.  In  the  following  we  explain  it  briefly.  At  the 

beginning  all  data,  including  those  referred  to  the  previous  time  step,  are 

passed  on  from  Box  I  to  Box  II.  In  Box  II  the  equations  for  the  vegetation 

layer  are  solved;  in  the  first  iteration,  as  a  zero  approximation,  the  values 

of  the  soil  moisture  and  of  the  temperature  at  the  surface  are  those  calculated 

in  the  previous  time  step.  Then  the  values  of  T,  ,  E,  ,  Ef  and  are 

f  ^  '  -tj*  9 

passed  on  to  Box  III  and  /  or  to  Box  I'J.  In  Box  III  the  equation  for  the  heat 
transfer  is  solved  and  a  temperature  profile  is  obtained.  This  is  passed  to  Box 
lU.  Successively  in  Box  IV  the  equation  for  the  soil  moisture  is  solved  and  a 
moisture  profile  is  obtained.  Finally  the  values  of  the  surface  temperature  and 
moisture  calculated,  r espec t i ve 1 1 y ,  in  Box  III  and  Box  IV  are  passed  on  to  Box 
II  and  a  new  iteration  starts.  At  the  end  of  the  iterative  process  a  new  time 


step  in  considered. 

In  Fig.  2  a  scheme  o-f  the  program  which  solves  the  vegetation  layer 
equations  is  given.  For  the  calculation  of  the  transpiration  rate  and  of  the 
evaporation  rate  and  the  sensible  heat  flux  at  the  ground  surface  some 
parameters  have  to  be  specified.  They  are  listed  in  APPENDIX  D  where  the 
listing  of  the  program  along  with  a  table  of  the  variables  are  presented. 

In  Fig.  3  a  scheme  of  the  program  which  solves  the  soil  moisture  floi/i 
equation  without  considering  the  effect  of  temperature  gradients  is  given.  For 
the  calculation  of  the  change  with  time  of  the  soil  moisture  distribution  an 
initial  profile  along  with  two  boundary  conditionsat  each  time  step  have  to  be 
specified.  The  boundary  condition  can  be  either  the  value  of  the  soil  moisture, 
pri»55ure  or  the  value  of  the  moisture  flux.  In  the  latter  case  at  the  first 
step  the  value  of  the  soil  moisture  presure  at  the  boundaries  his  to  be  chosen 
in  order  to  match  the  impposed  flux.  In  APPENDIX  E  the  listing  of  the  program 
along  with  a  table  of  the  variables  are  presented. 

For  the  tests  some  cases  for  wich  there  are  experimental  data  and  results 
of  calculations  carried  out  by  means  of  some  other  methods  have  been  picked  out. 
The  tests  confirm  that  all  Boxes  work  normally  and  that  the  soil  moisture 
computation  i saccompl i shed  successfully.  Some  results  of  our  soil  moisture 
computations  are  presented  in  Fig.  5  and  7.  They  show  the  change  with  time  of 
the  soil  moisture  distribution  with  depth  for  two  different  boundary  conditions 
and  for  the  same  values  of  the  physical  parameters. 

Details  on  the  physical  parameters  are  given  in  the  legends. 

Some  experimental  data  and  numerical  results  obtained  by  (1)  for  the  same 
case  considered  in  Fig.  5,  are  presented  in  Fig.  6. 

The  comparison  of  Fig.  5  and  Fig.  6  seems  quite  convincing.  At  last,  Fig. 8 


illustrates  how  Box  II  works.  It  shows  the  result  o-f  our  computation  ot  the 
diurnal  march  o-f  the  -foliage  temperature.  In  the  legend  o-f  Fig.  8  is  given 
whole  information  about  this  r^jn;  some  of  the  input  values  are  taken  from  [51. 

At  present  time  we  are  continuing  to  work  with  the  program.  We  are 
debugging  the  interaction  between  all  four  Boxes.  Simultaneously  we  are 
preparing  mass  tests  of  the  suggested  method. 


5.  COf^CLUSIONS  Ar4D  RECOMErJDATIONS 


In  conclusion  we  note  the  following.  A  method  which  allows  to  compute  the 
distribution  o-f  the  soil  moisture  and  temperature  with  respect  to  the  depth  at 
any  time  is  worked  out.  Also  the  moisture  and  heat  -fluxes  from  the  soil, 
covered  by  the  vegetation  into  the  atmosphere  can  be  computed. 

The  method  needs  the  following  information;  data  about  the  march  of 
meteorological  elements  at  the  shelter  level  (wind,  temperature,  humidity, 
precipitation),  data  about  the  march  of  the  cloudiness  and  solar  radiation,  data 
about  some  physical  parameters  of  the  soil  and  the  vegetation  cover  whose 
thickness  is  assumed  to  be  known  and  to  be  less  than  the  height  of  the  shelter 
1  eve  1 . 

The  method  is  based  on  the  numerical  solution  of  a  system  of  equations 
describing  the  heat  and  moisture  transfer  in  the  soil  and  between  the 
so i 1 -vege tat i on  system  and  the  atmosphere.  The  numerical  method  which  we  use 
consists  of  the  well  known  Crank-Nicol son  predictor-corrector  scheme  in  the  form 
suggested  by  Douglas  and  Jones  [191,  t203,  two  internal  and  one  external 
iteration  processes.  Such  a  method  provides  the  necessary  accuracy  and 
convergency  of  the  solution  and  it  is  rather  economic  from  the  point  of  view  of 
the  computer-'s  time.  At  present  time  the  programming  is  in  act.*^ 

The  theoretical  part  of  the  work  was  completed  and  accepted  for  publication 
in  "Progress  in  Desert  Research",  proceedings  of  the  Blaustein  Institute  for 
Desert  Research  (Israel)  -  UCLA  (U.S.A.)  Symposium,  1985. 

y It  should  be  noticed  that  although  the  formulated  problem  is  solved, 
however  the  work  is  not  finished  completely.  It  can  be  explained  by  the 
following,  for  some  reasons  which  are  beyond  the  author  control  he  was  almost 
deprived  of  a  programmer  or  of  any  assistance  during  his  work  with  the  contract. 
Only  in  the  last  few  months  he  had  a  programmer  for  one  day  work  a  week.  At 
last,  in  December  1985  the  author  has  received  an  assistance  by  a  qualified 
scientst  and  programmer  (M.G.  Scarpino)  who  had  an  opportunity  to  make 
^quaintance  with  the  problem  during  two  months  autumm  1985.  At  present  time  the 
work  is  being  carried  out  quickly  and  it  might  be  completed  in  a  few  months. 


Flowchart  tor  soil  moisture  and 
temperature  computation. 
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Program  Soil  scheme.  The  program  can  run  with  two  types  o-f  boundary  condition 
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A  scheme  the  calculation  o-f  sjitis-fyjng  to  the  boundary  condition  at  the 


surface. 
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The  distribution  ot  ^  with  respect  to  the  depth  of  the  soil  tor  t=0.1h,  0.2h, 
..0.8h  computed  by  our  program  tor  the  problem:  U.-p  =13.69cm./h  at  2=0,  •9' =0.1  at 
2= /0cm,  and  at  t=0 .  ~^Sec  ),  it  is  assumed 

d,  ){i  i-  0.  ^ 

where  k  -  34cm/h,  =0,287,  ^=0.075.  Field  temperature  and  moisture  remoual 

by  the  roots  ettects  are  neglected. 


0.2 
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FIG.  6 


The  result;  of  computation;  carried  out  by  the  different  numerical  methods  for 
the  same  problem  and  for  the  same  values  of  the  parameters  as  at  Fig.  5.  Also 


the  experimental  data  are  presented.  This  picture  was  taken  from  [1]. 


Calculated  diurnal  march  o-f  and  given  diurnal  marches  o-f  ,T^  and  F 
was  assumed  that  q  =0,007,  q  =0 .002+0 .005(  tK'6)  ,  u  =1  m/'sec,  C,.=0.05^ 
Cq  =0  .096 ,  ^.=0 . 2,  ^^j=0 , 1  .  All  this  information  was  taken  from  t53. 
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APPENDIX  A. 


DETERfllNATlON  OF  (x^).^  AND  UHETHER  CONDENSATION  OR  E'v>APORATI ON  IS 


OCCURRING  AT  THE  FOLIAGE  AT  t  =  . 


Setting  »*=  1  in  <3<5),  we  solve  the  eq.  F(X,  1)  =  0  numerically  by  the 
Newton-Raphson  method  (see  tl71): 


(55> 


The  iteration  process  is  stopped  when 


(54) 


tr 


The  tact  that  and  do  not  change  the i 

.fo) 

signs  at  _^  >  0  provides  the  convergence  o-f  the  iterations  -for  any  X  >0. 

At  ter  determining  )(  ,  we  compute  by  (37). 

Further  computation  depend  on  the  values  ot 

It  0,  then  in  accordance  with  (42)  condensation  on  the  toliage  is 

occurring  and  T-X. 


It  0,  where  ^  is  very  small  (we  put  ^  =  10  >,  evaporation 


trom  the  toliage  is  occurring.  However,  in  this  case,  the  evaporation  is 


su 1 1 i c i en 1 1 y  smal 1  to  neglect  it  by  putting  =  0 .  Then  we  t i nd  trom  (43) 
and  then  determine  trom  (14)  and  (17).  (Owing  to  (2),  is  negl  i  g  i  b  1  > 


smal 1  as  we  1 1  as 


s  .)  In  this  case  is  so  close  \.o  X  ,  that  it  is  possible 

It  '£[X)  >  ,  then  evaporation  trom  the  toliage  is  occurring.  IJe  denote 


to  take  ~T^  n  X  • 


One  can  show  that  then 


solution  of  the  eq, 


C/  f 


T'<  1  <  T' 


wher  e 


7-W  . 

//•  I 


0.  We  turn  again  to  the  Newton-Raphson 


method,  finding  T'^^Aron,  the  same  iterative  formula  (55)  in  which  one  should 
substitute  instead  of  1,  as  a  second  argument  of  the  function/^. 

But  now  it  is  better  to  set  ~  . 

A  criterion  for  stopping  the  iterations  is  the  same  as  (56). 

Proceeding  to  the  solution  of  the  eq.  (47),  we  note  that  the  function 

7-  .  -rfv 

is  continuous  and  increasing  at  the  interval  j  ^  5  ' 

of  (21),  (30),  (44)  and  (48)  one  can  show  that 


With  the  help 


j(T^]  -  1-8  I 

Consequently,  since  //*77)is  a  monotonic  function,  the  solution  of  the  eq. 
(47),  e.g.  If  //  /=  0 ,  exists  and  is  unique.  To  find  this  solution,  we  use  the 
so-called  "rule  of  false  position  or  reguli  falsi"  [19],  which  is  presented  b;- 


the  following  convergent  iteration  process 


y  " 


J - ;; - 

^  ^  rn; 


ij  j  <0 

u }  >0 


..h.r,  d,n,„  f'’. 


The  criterion  of  stopping  the  iterations  is  the  same  as  (56) 
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APPENDIX  B 


NUMERICAL  SCHEME  FOR  SOL'JIMG  THE  TEMPERATURE  EQUATION  (50). 


Since  C  and  A  do  not  depend  on  r,  the  eq.  (50)  is  linear.  The  boundary 
conditions  (27)  and  (53)  are  also  linear.  There-fore,  one  can  use  the  well  known 
Cr ank-Ni col  son  scheme,  which  is  absolutely  convergent,  to  provide  an  accuracy  ot 
the  order  o-f  and  is  amenable  -for  solution  via  the  so-called 

tridiagonal  algorithm  [1?]. 

Supplementing  the  time-grid  introduced  earlier,  we  shall  introduce  a 
uniform  space-grid  y  C=  0,1,... A/  and  denote  » 

luhere  ^  is  any  variable  ot  the  problem. 

Presenting  the  eq,  (50)  in  the  -form 

SSTL  4-  ZZ1A..2Z 

A  '~d-i  " 

we  replace  it  by  the  t  i  n  i  te-di -f  terence  equation  in  accordance  with  the 
Cr ank-N i col  son  scheme  and  then  transform  it  to  the  following  form: 


(60) 


(61) 


where 


ff 


(62) 


Eq .  <61)  is  solved  by  tridiagonal  algorithm  <5ee  tl*?]): 


=<>i  ~ 


<6  3> 


o(.  -TT  {0 


To  employ  this  aloorithm  one  needs  to  know  oC  i 

^  U 


and  . 


The  latter  value  is  known  from  <27)"7'  *7”. 

„/  /> 

To  find  Oc^  and  p  ,  we  have  to  rewrite  (53)  in  the  finite-difference  form. 

2^ 

Remember  that  the  Crank-Ni col  son  scheme  provides  accuracy  of  the  order  of 


with  respect  to  ^  .  Therefore,  we  use  the  following  representations 


(6*1) 


providing  the  same  accuracy. 

Substituting  i6A)  into  (53),  we  come  to  a  certain  linear  equation  linking 
V'  ,  *77  and  *77  .  The  second  linear  equation  linking  these  three  values  is  (61) 

0  T.  Ji 


at  i=  1.  Elimination  of  ~T1  from  these  two  equations  yields 


wher  e  o/,  •  and  /5 ^  are  values  of  interest: 


<.  65  > 


Mote,  that  in  the  expressions  <54)  tor  and  A”  the  sur-face  temperature"?^ 


/  - -  - - 'J 

reters  to  the  preuious  iteration  step.  The  re-fore,  the  computation  by  U'3) 


’'J 


can  be  considered  as  a  realization  o-f  one  iteration  step,  e.g.  as  obtainingY" 


Substituting  the  ualue”77.  obtained  from  <<i3),  into  (54),  we  proceed  to  the 
next  iteration  step,  etc.  As  a  criterion  -for  stopping  the  iterations  we  use  the 
i  n  e  q  u  a  1  i  t  y 


(67) 


APPEtJDIX  C 


NUMERICAL  SCHEME  FOR  SOLVING  THE  SOIL-MOISTURE  EQUATION  51 


First,  instead  o-f  we  introduce  a  new  independent  «Mr  table  ^"KirThhcH 
tr  ans-f  ormat  i  on" ) . 


The  eq .  (51)  takes  the  form: 


^ _ -ir  L. 

C-V^-r 


C-XKK  * 


The  terms  in  brackets  are  known,  because  Tf  or  the  moment j  has  been  -found  and 
in  the  expression  ‘6'  for^  soil  moisture  P  is  taken  from  thej-l  moment. 

Eq.  (69)  is  more  con'iement  for  the  numerical  solution  than  the  eq  .  (5r'> 
inasmuch  as  powers  in  coefficients  are  not  so  great,  although  the  parameter  ^ 
can  be  of  order  of  10  (see  Table  2>  from  tl2]. 

In  addition,  the  second  order  derivative  term  is  transformed  into  its 
simplest  form.  As  was  stated  earlier,  we  use  the  Douglas  and  Jones 
pr edi c tor -corr ec tor  method  1201  ‘or  the  solution  of  eq.  (69). 

First  of  all,  we  replace  (69)  with  the  system  of  two  finite-difference 


equations  (the  first  of  which  is  the  predictor,  and  the  second  is  the  corrector 
following  tl9]  (eqs.  (3-50a'  (3-50b)  in  1191).  Then  both  of  these  eqs.  are 
reduced  to  the  standard  tridiagonal  form. 


The  predictor  1=  gii.'en  b> 


<  i  2  ^  %  f 


-  F 

i 


where 


A.-i 


B,.  =  i 


A 


F=5A,,  u:  .  *^a..  -{v.  -V.  )y.f 
‘  v'  w-‘  ^di.rn  '-ij-i  4 


V 


/.  1-2  q  .-  4*--? 
i  —  _ ■{/ 


1  -  c-wi^  - 


ij  ii 

fX-a-f  ^Tj-i  ^  -i  ^ 


(70 


(71 


and  is  -followed  b>  the  corrector; 


y? 


•J 


where 


To  solve  (70)  and  (72)  we  use  the  same  -formula  (<43)  ot  the  tridiagonal 


algorithm.  It  should  be  noticed  that  eqs,  (70)  and  (72)  are  solved 
successively.  First  -from  predictor  (70)  we  -find  ^  ■  Substituting 

these  values  into  (73),  vje  solve  the  corrector  (72),  -findino  IX.  . 


The  values  and 


have  to  be  found  from  the  boundary  conditions  at 


the  surface  (7).  1*o  satisry  (7)  while  computing  XX-  ■  ,  we  emplo>  the  logical 

yj 

scheme,  presented  in  Fig.  4. 

It  should  be  not^ed  that  the  values  relate  formally  to  the 

intermediate  time  step  ,  but  they  can  be  interpreted  as  the  first 

i  2 

approximation  to  the  values  of  interest  XX"- .  .  In  this  regard,  the 

V 

coefficients  figuring  in  the  boundary  conditions  (7>  refer  to  the 
moment  of  time  in  the  predictor  as  well  as  in  the  corrector. 

The  bottom  boundar-.-  condition  (27)  gives 


V  ,~V  ^ 


^  (  1  for  dry  soil 


1  for  satur  ated  soi  1 


One  can  see  that  an  aoproximate  satisfaction  of  (7'  is  achieved  at  the  cost 
of  either  the  predictor,  or  the  corrector,  or  both  the  equations  hsving  to  be 
solved  twice.  At  worst,  the  computations  are  doubled. 

The  'lalue  U} .  obtained  from  the  corrector  is  taken  as  the  input  lalue  for 
the  next  i ter  at  ion. 

The  iterations  stop  when 


ynax  iX. . 

/ 


'i/7^  w  ’ 


This  criterion  provides  1";  accuracy  for  Q 


The  values  c/^and  corresponding  to  the  conditions  No  1  and  No  3  -folloi.i 
•from  the  last  expression  o-f  (<53)  at  L=  1)  rewritten  -for  .  lie  have  U.-O 

^=1.  To  obtain  these  parameters  -for  the  conditions  No  2  and  No  A  one 
must  satistx  the  moisture  balance  equation  at  the  surface  when  the  runof-f  is 
absen  t 


rewritten  in  the  -finite-difference  form  for  the  moments  "Zf •  Z  and  • 

J  S.  J 

IJe  carried  out  the  following  transformation  with  (.76). 

1.  Introduced  2/  instead  of  ^  by  (<48>. 

2.  Replaced  the  der  i ','at  i  ves  by  the  three-point  finite-differences,  like  (6A)  . 


After  that  the  procedure  of  obtaining  And  is  exactly  the  s 


ame  as  in 


tho  case  of  the  so i 1  -  temper  a tur e  equation. 


For  the  condition  No  2  we  have 


ex'  =  LJl!4iL 

I  *  nt 


Mi  ^  i-P. 

Q  =  —I - i-i- 

'‘t  I  -f-  m*  ^ 


where 


J  C-V)^r 

I  /p_£-  _  y) 


3.A2 


For  the  condition  No  <1  we  have 


The  express  i  on=  tor  the  runo-f-f  are  obtained  in  this  wa/: 


The  constant  dimensional  factor  in  -front  of  the  expression  for  R 
dropped.  One  must  determine  only  the  sign  of  this  expression. 
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APPENDIX  D. 


Tf^BLE  OF  THE  VARIABLES  OF  THE  PROGRAN  GREEN 


Physical  constants  and  so i  1 -'-'eoe tat  i on^atmosphere  data 


AD 

AHC 

LATC 

GRAV 

RW 

SIGMA 

MSAT 

HaULT 

FA 

FE 

GE 

3TRES 

FF 

NFflAX 

TRCB 

IJC  1 

UC2 

ANLE',» 

CAfJLEV 

SPNAX 


a  i  r  den  si ty 

speci-fic  heat  of  air  at  constant  pressure 

latent  heat  of  vaporization 

or av it"  acce 1 er  a t i on 

gas  constant  for  I'later  vapour 

Stef an-Bol tzman  constant 

moisture  value  at  saturation 

moisture  wilting  point 

f ol i age  albedo 

foliage  emissivity 

ground  emissivity 

generalized  stomatal  resistence 

area  averaged  foliage  shielding  factor 

maximum  liquid  water  retained  by  foliage  per  unit  horizontal 
ground  area 

dimensionless  heat  and/or  moisture  transfer  coefficient  for 
a  bore  soil 

i t  CO i nc i de s  w i th  u  (24)  when  T,  '  1- 

it  coi  nc  i  des  wi  th  u^  (24)  when  T^^  )  T^^ 
anemometer  height 

canopy  height 

maximum  downward  shortwave  radiation  at  height  "h" 


V' 


V 

V 
S 


I 

r. 

I 


ft  h 


S 


k 

i: 

« 


f* 


Physic  3l 


AD 

AHC 

LATC 

GPA',.' 

RUI 

SIGMA 

MSAT 

rWILT 

FA 

FE 

GE 

3TRES 

FF 

MFHA( 

TPCB 


UCl 

A 

UC2 

.> 

> 

A 

aMLE',' 

CAMLEV 

SPMAX 
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APPENDIX  D. 

TABLE  OF  THE  VARIABLES  OF  THE  PROGRAM  GREEN 

constant;  and  soi 1 -veQ9 tat i on-atmosphere  data 

:  air  den  =  i  tv 

:  rpecitic  heat  o-f  air  at  constant  pressure 

:  latent  heat  ot  vaporization 

:  gr av i ty  acce 1 er a t i on 

:  gas  constant  -for  water  vapour 

:  S t e +' an -Bo  1  t zman  constant 
:  moisture  value  at  saturation 
!  moisture  wilting  point 
!  -f  ol  i  age  al  bedo 
:  -f  ol  i  age  em  i  ;s  i  v  i  t  y 
!  ground  emissivitv 
:  generalized  stomatal  resistence 
:  area  averaged  -foliage  shielding  -factor 

:  maximum  liquid  water  retained  by  foliage  per  unit  horizontal 
ground  area 

:  dimensionless  heat  and/or  moisture  transfer  coefficient  for 
a  bore  soil 

:  it  coincides  with  u  (24)  when  T.  <  T- 

C  h  ^ 

:  it  coincide;  with  u,  (26)  when  T.  >  T. 

•  K  ft. 

:  anemometer  height 

:  canop-y  height 

:  maximum  downward  shortwave  radiation  at  height  "h" 


3 


'■'ariable;  appear  inq  as  input  '.'alues  at  a  Qi'.’gn  time; 


CF 

PR 

TA 

TH 

TF 


TG 


QA 


QG 


UA 

SR  HO 


WFII 


riPfllN 


MG 


:  cloijd  fraction 
:  pr ec  i  p i tat i on  rate 
:  air  temperature  outside  the  canopy 

!  air  temperature  at  canopy  height 

-foliage  temperature  (only  a  first  evaluation  to  begin  the 
compu  t  at i on) 

:  ground  temperature 

:  air  specific  humidity 
:  ground  specific  humidity 
:  mean  wind  up|ci..;t--  outside  the  canop>' 

!  downward  shortwave  radiation  at  height  "h" 

:  liquid  water  retained  by  foliage  per  unit  horizontal  ground 

area 

:  minimum  moisture  in  the  root  zone 
:  moisture  at  the  surface 


COTiput  at  i  on  variables! 


DTA  :  time  interval  in  the  finite  difference  equation  of 

conservation  of  the  viater  on  foliage 


Ml  Tri'( 


ma;rimum  number  of  iterations 


CF 


C 1 oud  f r  5r  t i on 


PR 

TA 

TH 

TF 


TG 


QG 


SR  HD 
WFM 


riRtilM 


tlG 


:  precipitation  rate 
:  air  temperature  outside  the  canopy 

:  air  temperature  at  canopy  height 

-foliage  temperature  (only  a  first  evaluation  to  begin  the 
coFTiDU  t  a  t  i  on) 

:  ground  temper at>  re 
:  air  specific  humidity 
:  ground  specific  humidity 

!  mean  i.jind  veiouit:--  outside  the  canopy 
:  downward  shortwave  radiation  at  height  "h" 

:  liquid  water  retained  by  foliage  per  unit  horizontal  ground 

area 

:  minimum  moisture  in  the  root  zone 
:  moisture  at  the  surface 


CofTiput  at  i  on  variables; 


DTA  :  time  interval  in  the  finite  difference  equation  of 

conservation  of  the  water  on  foliage 


MI  TMX 


maximum  number  of  iterations 
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Func  t i on  5 

;  "errion  of  the  Cl  3,u=  i  us-CI  apeyr  on  equation 

:  T:  temperature;  MG:  surface  moisture;  it  gi'.'es  the  '.'alue  of 
the  specific  humidity  at  the  surface  (25) 
as  defined  in  (37) 

Auxiliary  uariables 

R  :  fraction  of  potential  evaporation  rate  from  foliage 

PI  :  unitary  i,'?lue  of  P 

P2  :  '<a1ue  of  R  when  no  water  is  retained  by  the  foliage 
El  !  i  t  c 0  i  n c  i  de $  w i  t h  1  ,  ( 1  1 ) 

E2  :  i t  c  0 i n  c i de  5  w i t  h  f  2  (11) 

t 

E 3  :  i  t  c 0  i  n  c  i  de  5  w  i  t  h  S  3  ,  ( 1  2) 

AUXD  !  computational  constant 

PL  :  downward  longwave  radiation  at  height  “h" 

IjA  :  ground  albedo 

UC  :  wind  speed  it  coincides  with  M  ,  (.26) 

RIB  :  bulk  Richardson  number 

TPCH  :  dimensionless  heat  and/or  moisture  transfer  coefficient  at  the  top  of 
a  dense  canop> 

LIAF  :  mean  wind  velocity  within  a  canopy 
CFO  :  computational  variable 
CU  :  i t  CO i nc i des  w i th  C„  /  1 4) 

AU'!A  ;  computational  variable 


F03AT  (T) 

FOG  a,  110) 

FCECK  (T  )  : 


51 


AA  : 

it  coincid“;  i."ith  .''33> 

AF  : 

it  CO i nc i de =  w i th  a,  <33) 

AG  : 

it  coincides  with  a.  ,(33> 

AE  : 

0 

it  coincides  with  a^  ,<35) 

PE  : 

it  coincides  with  ^(38) 

'(lAUX  : 

computational  variable 

XI  ! 

it  coincides  with 

AUXF  : 

computational  variable 

AO  : 

it  coincides  with  A  ,  (48) 

BO  : 

it  coincides  with  B  ,''.48) 

OAG  : 

it  coincides  with  q  <35> 

TAG  ! 

it  coincides  '.•) i  t h  ,  ( 39 ) 

NITER  : 

number  of  iterations 

CECK  ! 

value  oi  ) ,  (37) 

FAUX  : 

computational  variable 

FI  ! 

computational  variable 

F2  : 

computational  variable 

FDER  ! 

derivative  ot  the  -function  ,r),(36)  with  respect  to  the 

temperature 

TFP  ! 

foliage  temperature  value  referred  to  the  previous 

i ter  at i xe 

TFl  : 

foliage  temperature  calculated  by  solving  equation 

( 36)  with 

TF2  : 

foliage  temperature  calculated  by  solving  equation 

<no  water  is  retained  by  the  foliage) 

(36)  wi th 

EF  : 

evaporation  rate  from  foliage;  it  coincidence  with 

in  (40 

FTFl  ! 

value  of  j 5  defined  <57)  with  T=TF1 

FTF2  ! 

value  of  (  ^  defined  (57)  with  T=TF2 

tep . 
=  1 . 
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QSAT  s  saturation  spec i tic  humidity 

OAF  :  mean  specific  humidity  inside  the  canopy^  (31) 

TAF  I  mean  temperature  inside  the  canopy,  (32) 

RTF  :  fraction  of  the  potential  euapotr ansp i r at i on  rate  as 
T^,(45> 

F3  :  ualue  of  the  function  F<Tj  ,r),(36> 

FTF  :  ualue-of  the  function  ),(5?) 


3  function  of 


PPQGPAM  GREEH 


■  IHPlJT -OUTPUT  • 


THIG  PROGRAM  SOLVES  THE  EHERGY  BALANCE  EQUATION  FOP  A  '.'EGETATION 
LAYER  AS  PARAMETERIZED  BY  DEARDOPFF 

UNITS  APE*  CM  MIN  CAL  Gl'AMS 


BEAL  LATG  r  MG  - MMILT ■ MRMIH • MSAT 

DATA  GRA'v'/?31  ..'-Ru  -  96l  .  SES/ •  LATC/585  ,  /  •  AHC  .Z-l  ■•AD/l  .26E-3  '• 
t SIGH A  /  .  S26E-10/ -FA/ . 2/ -  FEZ  - - GE/ .S5/-STPES/2 , / -UPMAX/ • 1/ » 
i TRCB  /  .  0057  ■  -  R  1  '  1  .  -  MUIL  T  *  1 

FUNG T IONS *x«x****«* ******************************* ************ 


FQSAT  '  X  )  =  .  S8E-Z*E:<P  '  1~  .  269*  <X-27S  ,16)  '  <  X-S5 . 36  )  > 

PQG  (  Y  .  Y  ■>  =FGSAT  <  X  )  *EXP  <  GRAU*Y/  <  RU*X  '■  i 

FCEGK’'X)  =  AHC/LATC*'TAG-X)  +  1/<LATC*AE'*'RE-E3*X**^> 

x*:«x  *******************************  *************x**x*x******** 

SET  VEGETATION  PARAMETERS  AND  COMPUTATION  STEPS 

PBtht* -' ASSIGN  FF  ANLEV  CANLEV 
^'E  AD*  -  FF  .  ANLEV  -  CANLEV 
DTA  =  S 
NITMX=20 

AUXILIARY  CONSTANTS 


El=FF*GE  rSIGMA/ ^FE  +  GE-FE*GE) 

E2=';  1  -FF)*SIGMA*GE  +  E1*FF 

E3=E1+SIGMA*FE 

AUXD= 1’. 269*23^. 3*. 33E-2 

INPUT  VALUES  FOR  THE  VEGETATION  LAYER  EQUATION  SOLUTION 

PRINT  If  •' ASSIGN  CF  PR 
READ* -CF -PP 

PPrr'T*-  '  ASSIGN  TA  TH  TC  ' 

PPAP*  -  TA  -  TH  ■  rr- 
PPINT* •' ASSIGN  ha  nn' 

PFADt •QA-QC 
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PPTHT* ASSIGN  UA ' 

READ*  *11 A 
DC  I  =  SCO 
!irC=-6(f.0C 

PRINT  ASSIGN  SPMAX  SPHD  ’ 

READ* • SRMAV . SPHD 
PRINT* ASSIGN  WPN' 

READ* • WFN 

PRINT*.  ' ASSIGN  HPMIN  HG  HSAT ’ 

READ* .MPHIN-HG-MSAT 
TF  =  TA 

PL=  '  CF  +  1 . 2  1  *  <’  1  -CF  >  *QA**  *  08  >  *SIGMA*TA**4 
IF  <  HG  .  LE  .  HSAT  '■  GA=  ,  31  -  .  1  7*HG/MSAT 
TF''HG.GT.MSAT>  GA=,14 
IFdH.LT.TA)  LIG  =  LIC1 
IF^TH.GE.TA)  !JG  =  UC2 

PIE:  =  GRAU*  (  ANLE',>-FF*CANLEU  •>  *  <  1  -TH/TA  1  /  ( IJA**2  +  IJC**2  ) 

IFfRID.LT.O.  >  TRGH  =  AUXl*(l+24,f.*f-AUXl*RIE:^**.?' 

IF  fRIB.GE.O.)  TRCH  =  AUXl/(l  +  ll  .5*.PIE:> 

UAF=  (  1-FF+  ,  83*FF*TRCH**  ,5''».Uh 
CFrj=  .  01*  f  1  +  1800 /UAF 
CU=7.7E-2*(UAF+30^ 

TPCG='’  1-FF’>*TPCE:  +  FF*TPCH 
AUXA  =  TRC.G*.IIAF+FF*(CU  +  TRCH*IIA' 

AA=TRCH*UA/AUXA 
AF  =  FF*CLI  '■'ALIXA 
AG=1-AA-AF 
AE  =  AD*CLI*  >'  1  -AF  ' 

PE=  f l-FA)*SPHD  +  FE*PL  +  El*TG*»4 

XIAUX=1 / ( SRHAX  /  ( SPHD* . 3*SRHAX  ^  +  ( MUILT  /MRMIN  >  **2  ' 

XI  =  YIAUX  /  <  STPES*CFO*IJAF  ' 

R2  =  XI/n+Xi:' 

An  =  DTA*AD*CU*FF*  ^ 1 -FA  +  XI ) /UFHAX 
E:ri-UFN  +  DTA*FF*PP  'NFMA'’ 

0AG=<'  AA*QA  +  AG*nG'>  '  ■'  l-AF  ' 

TAr-=f  AA*TA  +  AG*TG  >  '  ■  2-AF  ' 

PRINT*. 'OAG  TAG  '.QAG.TAC 
AL!XF=-XI  1  -AF  +  XI^ 
r 

r  iiEMTON-PAPHSON  METHOD  FOR  THE  SOLUTION  OF  THE  CANOPY  ENERGY  E'.ALANGE 

C  EPUATION 

r 

C  ITERATIVE  CALCULATIONS 

C 

I  NITEP=NITEP+1 

IF  ^NTTER*EC'.NITHX'  GO  TP  " 

CECK  =  FCECKn  F  ' 
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FAUVrrFpqAT  '  TT  ’-PAr 
F  l  =  AUXD*FnSHT  '■  IF'  ■'  '  TF-  .  S6  ■) 

F2=  AHr  +  ^!*'F3»TF*»:3  'At  >  ‘'LATC 
FDEP  =  F1  +F:: 

TFP  =  TF' 

TF  =  TFP-.-FAUX-CECK)  '<F1+F2  ! 

TFC  APFaF-IFP' .GT.  .  01  •  GO  TP  1 
PPIHT»  ■ FIPST  COMPUTATION  TF='-TF 
TFl =TF 
C 

C  CECF  FOP  CONDENSATION  OP  NEGLIGIBLE  EVAPOPATION  CASE 

C 

IF- CECK.GE. .00001 >  GO  TO  2 
C 

C  CONDENSATION  CASE 

C 

PPINT*  .  ’  CONDENSAtlON  CASE*  TF=’  -TF 
FF  =  FF*AE)rCECF 

UFN  =  UFN+  ■'  FF!rPP  +  EF  )  »:DTA/MFMAX 
IF  (NFN.LE. 1 )  GO  TO  6 
PRINT* »'NFN  EXCEDED’ 

GO  TO  e 

r 

C  EVAPOPATION  CASE 

r 

2  P  =  P2 

NITEP=0 

TF1=TF 

?  NITEP=NTTEP+1 

IF  (NITEP.EQ -NITr;-' -  GO  TO  7  . 

CECf  =FCE  CF  '•  TF  ■ 

FAIJ'^  =  AljyF»  ^FOSAT  -'TF  --OAG) 

FI  =AUyp*AUXD»FQSAl  iTF>  '  f  TF-St- . ’=' ■  **2 
f:=  AHC  +  ‘!*E3*TF*»3  AE  ■  .-'LATC 
FDFP^Fl +F: 

TFP=TF 

TF  =  TFP-  '■  F  AU  y-CECI'  )  'FDEP 

IF  f  ABS  T F - TFP  r  -  C T  .  .  0 1  '  GO  TO  3 

PPINT* EVAPOPATION  CASE  FIPST  EVALUATION  OF  TF :  ' • TF 

tf:'  =  tf 

C 

C  METHOD  OF  PULE  OF  ' PEGULI  FALSI'  FOP  THE  SOLUTION  OF  THE  FOLIAGE 

r  MOISTURE  CONSERVATION  EOUATION 


7T  >■  EiO  .  LT  ,  ■  OOOOl  ■'  CO  TO  5 
NTTEF=0 

rTri.  =  1.-E:n  +  DTA»:AE»:rF*FCE:Cf' 'TFl  ;  'i-'FMAX 
FTF2=-B0 
TF=' TF1  +  TF2) 
fJITEP  =  NITEP+] 

IF  ■ NITER. EO.NITHX'  GO  TO  - 
CECF  =  FCECKnF  ' 
nSAT  =  FnSAT  TF  '< 

0AF=QAG+AF*CECF 
PTF=  < 1-AF  ^  *CECK/ < QSAT-QAF  > 

UFN=  !  (  XT  +  1  'j  »:PTF-XI  ^  .  C 

IF  '  MFN.LE.  1  .  )  GO  TO 
PRINTS- t'NFN  EXCEOEO' 

GO  TO  8 

F3  =  AUXF* ( OSAT-OAG  ^ -CECF 
FTF  =  UFN-A0*F3-E:0 
TF P=TF 

IF  ■'  FTF  .  GT  .  0  .  '  TF=  '■  TF2*FTF-TFR»:FTF2  ^  /  <  FTF-FTF2  ) 
IF'  FTF  ,LT  .  O'  TF=  f  TFP»FTFl-TFl5rFTF)  /(FTFl  -FTF  • 

IF  '  ARSdE-TFRl  ,GT.  .n  GO  TO  ^ 

PRINT*  .  •EVAPORATION  CASE  TF=’.TF 
T  AF  =  AA*Ti^i-tAF*TF'tAG*TG 
CiAF  =  OAG  +  AF»FCECr' TF  '. 

PRINT*'-'  TAF  OAF  '-TAF-OAF 
TH=  .1  -FF  *TG  +  FF  *TAF 

PRINT*  •'  '.'hLUES  HEHOPIZED  FOR  SUCCEBBI'.'E  BALANCE 
PRINT*  •  'TINE  TF  TH  NFN  '  -  TINE • TF f TH r NFN 

PRINT*- 'TIME  TG  TA  TF  ' - T IME • TG • T A » TF 
GO  TO  P 

PRINT*- ' NITER  EXCEDED  '-NITER 
CONTINLIP 
STOP 
ENO 
OF  file 
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APPENDIX  E. 


TABLES  OF  THE  VARIABLES  OF  THE  PROGRAM  SOIL 


Soil  d? 1 3 


Votum^tric  fiater  content  at  the  saturation, 
residua)  "olumetric  water  content 
initial  constant  no) ume tr i c  water  content 
saturated  soil  hydraulic  conductiuit> 
it  coincides  with  ^  in  the  -first  of  (12) 

it  coincides  i.»ith  A"^  in  the  first  of  (12) 

it  coincides  with  ji  in  the  second  of  (12) 

it  coincides  witlt  ^in  the  second  of  (12' 


Computation  uariables 


I  PR  If  IT 


number  of  layers 


number  of  levels 


depth  of  a  layer 
time  step 


2T/PT  i-ihere  T  is  the  duration  of  the  infiltration  process 
inde::  for  pr  inti  no 

computation  step:  there  are  two  computation  steps  for  each  time  step 


-  ^  ^ 

Variables  computed  at  each  iteration 


H':  1) 

soil  water  oressure 

TET<i:' 

'v'oi  ume  tr  i  c  water  content 

K(  i:' 

soil  hydraulic  conductivity 

C':r.> 

soil  water  capacity 

ck  I ) 

soil  water  diffusivity 

HOLD  ':i  :> 

soil  water  pressure  referred  to  the  peviou 

F'i  D 

element  ;f  the  t(  [diagonal  matrix 

0'  1' 

element  of  the  tri diagonal  matrix 

AC  ■;  1) 

• 

coefficient  of  the  tri diagonal  algorhythm 

EC-:  I ' 

« 

( 

coefficient  of  the  tri diagonal  algorhythm 

TIME 

. 

t  ime 

Anx  i  1  i  ary  ar  i  ables 


FACT 

:  computational 

aux i 1 i ar  / 

uar i abl e 

AU-'l 

:  computational 

au  ;<  i  1  i  ar  y 

ar  i  ab  1  e 

AUX2 

:  computational 

aux i 1  i  ar  y 

uar i abl e 

AU;''3 

:  computational 

aux i 1 i ary 

var i abl e 

F'jnc  t  i  ons 


FK<H' 
FC' H> 
FTET^H> 


soil  hydraulic  conductivity  as  a  -function  of  the  water  pressure 
soil  water  capacity  as  a  function  of  the  water  pressure 
soil  water  characteristic  function. 


-  •< 


-J" 


ft- 


'j 

.1 

t 

V 

k 

i 

i 

L- 


7 

i 

f, 


These  are  e^presions  ',12'  in  [1] 


O  1-^  o 


-r,9  - 


r- P 0 G p  AH  G  0 1 L  -  I i '! P U T  •  0 IJ T P U T  ■ 

C 

C  nHE-DIHEHSIOHAL  IHFILTPhTIOH  MODEL 

G 

r  Hh'>EPPAHP  PQIL  PIJNCTIOHS?  case;  sahd 

r 

C  LIMITS  ape:  cm  hour 

c 

r 

PEAL  KS 

DOUBLE  PPECISIOM  FP  •  PC  •  FTET  •  K  •  C  •  CK  •  F  .  G  .  AC  •  BC  »  H  r  HOLD  ?  AUX 1  •  ALIXZ- 
$ALIX:: 

DTMEMBTOhi  H  '  R  1  ■  .  P  <  F;  1  ■  •  C  '  B 1  )  .  CK  ''  B 1  )  .  G  ^  RO  '  •  F  ■'  RO  >  •  AC  ^  RO  )  .  PC  '  Rri  ,  . 

«  HOLD  '  fM  '  .  TFT  ^  P  ]  ■ 
r 

r 

DATA  TETR/0.0'5  •  TETS/0 . 287/ •  ALFI/0 . 621E-6/ ■  BETA ''3  .  Re.  -  f 
«  AT  ^0 .851E-6  /  *0  ''13;  69  •'  rDZ  M  ,  '  r  TETN  '0.1  /  .KS/39  .  •'  • 
tB/9 .7^.-  ,HH..  -61 .5 /’H  •'70/  .H1/-20.73/ 

C 

C  SOIL  F UN CT TO HS »■»:*»:** IK****** 

FK ( X  ^  =KS/ ( 1 + AI*ABS  <  X ; **B ' 

FC  f  V  ^  TETS-TETP  )  * ALFI*BET A*ABS  ( X  )  **  (  BETA-  1 )  / 

1  +ALFI*ABS(>';**BETA>**2 

FTET  ■'X:'  =  '' TETS-TETP  )  /a  +  ALFI*ABS(X)**&ETA)+TETP 

r  ******************************  X  *******  x;  *  X  ************************  * 


PPIin  *•  ’ ASSIGN  UMAX ' 

PEAD».JMAV 
DT=1 . /^20 
PM  =  N+  1 

FACT  =  D2**2.'DT 
r 

r  BOUNDARY  CONDITTON-  IMPOSED  FLUX*.  11  =  2:  FIXED  L'ALUE ;  IT  =  3 

F' 

IT  =  2 

r  11=:? 

r 

C  INITIAL  CONDITIONS 

C 

DO  10  1=1 -Ml 
H( I '=HN 
H0LD'1'=HN 
K' T''=FP':HN  ■ 

C  < 1 ■ =FC ' HN> 

TETf  T'-  =  FTET'HN  ' 

10  CPdl^Cd'/pfi; 

C - 


n 
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C  riPST  LE','F.L.  WHEi''  A  FLUX  IS  IHPOBED 

r 

PPIrn  »;■' ASSIGN  HI' 

REAF-it  -  HI 
H'  1  ''-HI 
HOLDC  ]  .  .=:HI 
K."  1  ^=FKfHI) 

C 

C  FIRST  LEVEL  FIXED 

r 

r  H':n=Hi 

r  HOLr'(i''=Hi 

r 

PRIHT* INITIAL  TETA  PROFILE' 

PRINT  100- Vl.TFT'T*  rl=l  .N1  ' 
r  PRINT*.  '  INITIAL  RRSSURE  HEAD* 

C  PRINT  100.  (I.Hd  ■  .1  =  1  rNl  ) 

PRINT*. ■ INITIAL  VALUES  K  C  CK ‘ 

PRINT  200-  ■' I  •  P  i  I  '.  C '■  T CK  d '.  T=1  .  N1  ■> 

r 

J  =  ■ 

IPRINT=1 

TIME=0. 

r 

C  EIEC-INNING  TIME  LOOP 

r 

1  .J=  J+i 

TF'-J  2*2.NE.J"'  GO  TO  2 
TIHE=TIHE+DT 
r 

T'  n  2  1>  2  ■  f  ' 

AUX  1  =  '  K  '  J-^l  '  d-1  '  'K  '  I  • 

AIJX2=  '  H  d+  d  -  H  d-  1  '  ^  /  2 

ALIX3  =  H0L  D  d-*  1  ■  -2  »  HOLD  d  ' +HOLD  d- 1  • 

I"'  F  '■  I  ■>  =2»CF' d  '  *FACT*HOLD  d:  +AIJX1  *  ^  AUX2-D2  >  +AU 

r 

!-,n  TO  2 

r 

2  DO  20  1  =  2..  N 

ALIX  1  =  (  K  dl-*  3  -t-  (  I-  3  •'  >  'K  ■'  I  . 

AIJX2=  '  M  d  -^  1  >  -H  d  -  1  ■>  '  '2 
30  F  d  ■'  =2*CK  d  •'  *H  d  ■  »FACT+-AIIX1  *  <  AUXE-DZ  *  '2 
C 

3  CONTINUE 
C 

DO  TO  1  =  2 -N 

TO  G  (  T  '  =2*  (  3  -^F  ACTxCK  '  I  '  • 


-  fll  - 


riXETD  '-'ALUE  AE  EOMMUAEY  roMDITIOM 

r 

r  hC  ■  IT'-  =  1  'G' 2  ' 

C  E:r  (IT  '■  r  HI  +F  '  2  -  G  '  2  ' 

r 

2  IHF'OSED  FLUX  AF  GGUNDAF'Y  COHDITIOH 

r 

AC^IT  .=2-G(2  '  '2 

GC  (  IT  ^=F  (  2  ^  '2  +  D2*:  Q  ''K(  1  • 

r 

L 

DO  GO  Ir:TT-^l.fn 

Ar(I^  =  l  '(G(I-1  ''-AC<I-1  ■>  •' 

50  GC  (I  >  =  '  E;C  ( I-  1  )  +F  '  I-  1  ■  )  '  G  ( I-  1  ■>  '  AC  ( I-  1  >  ) 
r 

2  Ir'  THE  case  OF  A  FIXED  OALUE 

C 

2  DO  60  I  =  lTrJ-l 

C 

2  IH  THE  CASE  OF  IHF’OSED  FLUX 

C 

DO  60  1=1. H 
13  =H-  T-‘  ] 

HOLDdl  >=Han 

H(I1  '=AC(I1  +  1  >»:H(I3+1  ■•+E;r(Il+l  ' 

T'  1 1  ''  =  F  F  I"  H  '  I  1 
C  T  .1  '  =  F  r  (  H  '■  1 1  ■  ■ 

TET'll  '=FTET''  H  a  1 
60  CL’  Il  ■'=0  .'13  ■  'Kdl  ' 

2  6  IF  (  J/2!»:2  .  HE  .  J  )  GO  TO  .66 

6  IF  (  0  ,  EQ  .  3  ^T»  IF’F;I3JT  '  GO  TO  65 
GO  TO  66 

2 

65  F'FiIHT*  •  '  TETA  F’FOFILE  TTME^'-TTME 
F’F  THT  3  00.  d  ■  TFT  d  ■  ■  T=  3  •  M3  ■' 

2  F-  F’  I M  T  '  F’  F'  E  S  S  U  F'  E  H  E  A  D  F’  R  0  FILE’ 

2  -  F’F’TMT  3  00  '  d  •  H  d  •  d=  3  .  M  1  ■ 

TF'F'IriT  =  IF’F’I^iT-‘ 1 

1  RF’TMT».'3  C  Cl' 

2  F'  F’  TT!  T  ZOO.  d  •  F  d  ■  2  d  '  •  2 1'  dX'  •  I  ~  I  •  H 1  * 

66  IF  •' J  .  EO  .  UMAX  >  GO  TO  -q 

GO  TO  3 
T  20HTIMUE 

3  00  FOPHAT  >  5  '  T)  dl  .  E  3  3  ■  ■T  '■  ■ 

200  FOPHATf  ',I2.3E11 
STOP 

end 

--EOP-- 
END  OF  FILE 


